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ABSTRACT 



Context. Photospheric motions shuffle the footpoints of the strong axial magnetic field that threads coronal loops giving rise to 
turbulent nonlinear dynamics characterized by the continuous formation and dissipation of field-aligned current sheets where energy 
is deposited at small-scales and the heating occurs. Previous studies show that current sheets thickness is orders of magnitude smaller 
than current state of the art observational resolution (~ 700 km). 

Aims. In order to understand coronal heating and interpret correctly observations it is crucial to study the thermodynamics of such a 
system where energy is deposited at unresolved small-scales. 

Methods. Fully compressible three-dimensional magnetohydrodynamic simulations are carried out to understand the thermodynamics 
of coronal heating in the magnetically confined solar corona. 

Results. We show that temperature is highly structured at scales below observational resolution and nonhomogeneously distributed 
so that only a fraction of the coronal mass and volume gets heated at each time. 

Conclusions. This is a multi-thermal system where hotter and cooler plasma strands are found one next to the other also at sub- 
resolution scales and exhibit a temporal dynamics. 

Key words. MHD — Sun: corona — Sun: magnetic topology — turbulence — compressibility 



1. Introduction 

Magnetohydrodynamic (MHD) turbulence has long been impli- 
cated as a key process in coronal heating (Eina udi et al.L fl996) 
as well as fast and slow solar wind acceleration. Numerical sim- 
ulations have been crucial in investigating the role of MHD tur- 
bulence. 

To tackle the turbulence problem correctly with direct nu- 
merical simulations, it is necessary to resolve enough spatial 
scales to obtain an energy containing range, an inertial range, 
and a dissipation range. Because of storage needs, the more the 
complexity of the problem is reduced, the more spatial scales are 
obtainable. Hence most previous research in this area has been 
restricted to d ynamical effects, e.g., reduced MHD (RMHD) 
(IStraussL Il976l) and cold plasma models. This reduces the size 
of the problem considerably. For example, for RMHD the dy- 
namics of the system are determined by two var iables: a stream 
funct ion and a magnetic potential function dRappazzo et all 
I20071) . For the cold plasma model it is customary to evolve 
three velocity field components and three magnetic field compo- 
nents, neglecting the pressure term (thus not correcting for the 
solenoidality of the ve locity field) so that six variables are used 
dDahlburg et al.Ll2009t) . 

In all these approximations all thermodynamic al aspects are 
lost, i.e., there is no information on the connection between the 
dynamics of the system, responsible for the heating mechanism, 
and the temperature and density behavior. In particular energy 
lost through Ohmic and viscous diffusion is simply lost from the 



system and the thermodynamic physics involved with thermal 
conduction and optically thin radiation is neglected. 

This lack of information prevents a complete reproduction of 
the energy cycle: photospheric motions shuffle magnetic field- 
lines footpoints injecting energy into the Corona and trigger- 
ing non-linear dynamics tha t form small-scales, organized in 
current sheets (lParkerlll972h ,where energy is then transformed 
into thermal (and kinetic and perturbed magnetic) energy by 
means of magnetic reconnection. Heat is then conducted from 
the high temperature corona back toward the low temperature 
photosphere where it is lost via optically thin radiation, but at 
the same time causing chromospheric evaporation. 

Thermodyn amics are typ ically studied with ID hydrody- 
namic models dRealel [2010) that study either the equilibrium 
condition of a heated loop, or the time-dependent evaporative 
and possible condensation flows. However they introduce an ad 
hoc heating function to mimic the energy deposition and neglect 
the cross-field dynamics (including field-lines reconnection). It 
is therefore crucial to use the computationally more demanding 
three-dimensional compressible MHD (with eight variables), in- 
cluding an energy equation with thermal conduction and the en- 
ergy sink provided by optically thin radiation. 

RMHD turbulence simulations dRappazzo et all 120071 120081) 
have shown that in the planes orthogonal to the DC magnetic 
field the length of current sheets is of the order of the pho- 
tospheric convective cells scale (i.e., ~ 1,000 km), while their 
width is much smaller and limited only by numerical resolution. 
In order to understand the thermodynamics of such system in 
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this letter we carry out simulations that resolve the scales below 
the convection scale. The boundary conditions used in this pa- 
per do not allow evaporative flows to develop along the heated 
loops, our focus here is on the relationship of turbulent current 
sheet dissipation and the forming coronal temperature structure. 

Our new compressible code HYPERION has allowed us 
to make a start at examining the fully compressible, three- 
dimensional Parker coronal heating model. HYPERION is a par- 
allelized Fourier collocation finite difference code with third- 
order Runge-Kutta time discretization that solves the compress- 
ible MHD equations with vertical thermal conduction and radi- 
ation included. The results provided by this new compressible 
code preserve the spatial intermittency of the kinetic and mag- 
netic energies seen in earlier RMHD simulations, but also pro- 
vide new information about the evolution of the internal energy, 

1. e., we can now determine how the coronal plasma heats up and 
radiates energy as a consequence of photos pheric convection of 
magnetic footpoints (Dahlbur g et al.Ll2010l) . 

In Section 2 we describe the governing equations, the 
boundary and initial conditions, and the numerical method. In 
Section 3 we detail our numerical results. Section 4 contains a 
discussion of these results and some ideas about future directions 
of this research. 

2. Formulation of the problem 

2.1. Governing equations 

We model the solar corona as a compressible, dissipative mag- 
netofluid. The equations which govern such a system, written 
here in dimensionless form, are: 
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with the solenoidality condition V ■ B = 0. The system is closed 
by the equation of state 



p — nT. 



(5) 



The non-dimensional variables are defined in the following 
way: n(x, f) is the numerical density, v(x, f) = (u, v, w) is the flow 
velocity, p(x, f) is the thermal pressure, B(x, f) = (B x , B y , B z ) is 
the magnetic induction field, J = V x B is the electric current 
density, 7\x, f) is the plasma temperature, = p(djVj + d(Vj) - 
AV ■ \6ij is the viscous stress tensor, <?y = (dj\>i + d(Vj) is the 
strain tensor, and y is the adiabatic ratio. To render the equa- 
tions dimensionless we have used «o, the photospheric numer- 
ical density, Va, the vertical Alfven speed at the photospheric 
boundaries, Lq, the orthogonal box length (= L x - L y ) and To, 
the photospheric temperature. Therefore time (f) is measured in 
units of Alfven cross time (ja = Lq/Va). 

The function A(T) that describes th e temper a ture d epen- 
dence of the radiation is taken following Hildnei] d 19741) . nor- 
malized with its value at the base temperature To = 10000 K. 



The term Cm denotes a Newton thermal term which is enforced 
at low temperatures to take care of the delicate chromospheri c 
and transition region energy balance dDorch & Nordlundl fl999h . 
We use C N = 0.1 [r,(z) - T(z)]e- 2(z+0 - 5 « at the lower wall and 
C N = 0.1 [Ti(z) - T(z)]e- 2(0 5 L ~- z> at the upper wall, where Tfc) 
is the initial temperature profile. 

The important dimensionless numbers are: S v = 
PqVaLq/p = viscous Lundquist number, S = poVaLo/tj 
I . unciqui st number, B = popo/Bi = pressure ratio at the wall, 

5 /2 

Pr = C v p/kT ' = Prandtl number, and P ra d, the radiative 
Prandtl number plT 2 A n^K{To), determines the strength of the 
radiation. C v is the specific heat at constant volume. The term 
'V denotes the thermal conductivity. The magnetic resistivity 
(77) and shear viscosity (p) are assumed to be constant and uni- 
form, and Stokes relationship is assumed so the bulk viscosity 
^ = (2/3)//. 



2.2. Initial and Boundary Conditions 

We solve the governing equations in a Cartesian box of dimen- 
sions 1 x 1 x L z , where L z is the loop aspect ratio (0 < x,y < 
1, -L z /2 < z < L z /2), threaded by a DC magnetic field Bo in 
the z-direction. The system has periodic boundary conditions in 
x and y, and line-tied boundary conditions at the top and bot- 
tom z-plates where the vertical flows vanish (v z = 0) and the 
forc ing velocity in the x-y plane is impo sed as in previous stud- 
ies (iHendrix et all Il996t lEinaudi et all Il996t lEinaudi & Vellil 
I1999T evolving the stream function 



n I nt\ 1 nt n\ 

if>(x,y, t) = /, sin 2 _ + f 2 sin 2 ^— + - J , (6) 

where /,(x,y) = a' nm sin{k n x + k m y + C nm ), and v = x e ; . 
All wave-numbers with 3 < (A: 2 + A: 2 ) 1 ^ 2 < 4 are excited, so 
that the typical length-scale of the eddies is — 1 /4. As the typ- 
ical convective cell size is ~ 1 , 000 km and we use L z — 5 this 
implies that our computational box in conventional units spans 
4, 000 2 x 20, 000 km 3 . Every f, the coefficients a' nm and the £ nm 
are randomly changed alternatively for eddies 1 and 2. The mag- 
netic field is expressed as B = Boe z + b, where A is the vector po- 
tential associated with the fluctuating magnetic field b = V x A. 
At the top and bottom z-plates B z , n and T are kept constant at 
their initial values Z?o, no and To, while the magnetic vector po- 
tential is convected by the resulting flows. 

In what follows we have assumed the normalizing quantities 
to be: no = 10 17 »j- 3 , T = 10 4 K, Bo = 10~ 2 tesla, and L Q = 
4 x 10 6 m. It follows that the values of the various parameters 
appearing in the equations are: B = 1.7xl0 -4 , va = 6.9xl0 5 m/s, 
5.8 s, S = 2.7 x 10 9 , S v = 2.09 x 10 9 , Pr = 1.52 X 10" 2 , 



ta 



P rad = 3.35 x 10" b , k x r„ /z = 0.18 Wm-'K- 1 , In A = 10. 

The normalized time scale of the forcing, f , is set to 51.7 
to represent the 5 minutes typical photospheric convection time 
scale. The normalized photospheric velocity is Vo = 1.45 xl0~ 3 , 
corresponding to 10 3 m/s. 

We impose as initial conditions a temperature profile with a 
base temperature To = lO 4 ^ and a top temperature 8 x 10 s K 
with the following z-dependence for the nondimensional tem- 
perature T/(z) — 1 + 79 cos(nzlL z ), while the numerical density 
is given as n,(z) = l/r,(z). This choice allows us to neglect the 
gravitational effects since for most of the loop the temperature 
remains high enough to make the gravitational length-scale big- 
ger than L, at all times. 
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Temperature 
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Fig. 1. Temperature contours at t — 300 in the midplane (z = 
0). Maximum (130) and minimum (49) temperatures correspond 
to 1.3 million and 490 thousand degrees Kelvin in conventional 
units. Distances are normalized by Lq — 4 x 10 6 m. 



3. Results 

With previous definitions equation © can be replaced by the 
magnetic vector potential equation: 



dA 1 

— = vx(B e, + VxA) + -VxVxA 

at S 



(7) 



We solve numerically the equations (l)-(3) and (7) together 
with equation (5). Space is discretized in x a n d y w ith a 
Fourier collocation scheme (Dahlburg & PiconeL Il989h with 
isotropic truncation dealiasing. Spatial derivatives are calculated 
in Fourier space, and nonlinear product terms are advanced 
in con figuration space. A second-order centra l difference tech- 
nique (Dahlburg, Montgomery & Zang, 1986) is used for the 
discretization in z. A time-step splitting scheme is employed. 
The code has been parallelized using MPI. A domain decompo- 
sition is employed in which the computational box is sliced up 
into x—y planes along the z direction. 

We present the results obtained running the code with S = 
S v = 4 x 10 4 with a resolution of 128 3 . In order to keep the 
efficiency of the radiative and conductive terms in the energy 
equation as in the real corona, we have rescaled Pr and P ra d 
accordingly with the choice of S v , i.e., Pr = 793.39 and P ra d = 
0. 175. This choice is du e to the result found in the RMHD model 
( Rappazz o et al.L [2008b that turbulent dissipative processes are 
independent of viscosity and resistivity when an inertial range is 
well resolved. 

The system starts out in a ground state, threaded by a guide 
magnetic field: the interior of the channel has zero perturbed 
velocity and magnetic fields. The footpoints of the magnetic 
field are subjected to convection at the z boundaries with initial 
number density and temperature profiles assigned as described 
above. The behavior of the volume averaged quantities, such as 
kinetic and fluctuating magnetic energies and resistive and vis- 
cous dissipation show a tempora l behavior similar to the pre- 
vious RMHD results obtained bv lRappazzo et alj d2007l I2008L 
2010). 
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Fig. 2. Maximum temperature (T max ) and maximum electric cur- 
rent \j\ max as functions of time. 



The fluctuating magnetic and kinetic energies (e,. = |(«v 2 ) 
and eb = ^-(b 2 )), as shown in previous RMHD model, and the 
total internal energy (U = (nT)/3/(y - 1)), not obtainable in the 
RMHD case exhibit an intermittent behavior. After the system is 
driven by photospheric motions for more than 100 Alfven times 
a statistically stationary state is attained and all quantities fluctu- 
ate around their steady values. 

We have verified that the dynamics of this problem are well 
represented by RMHD. We define the volume-a veraged veloc- 
ity 6 V and magnetic 9b anisotropy angles as in IShebalin et al.l 
d!983l) . For isotropic turbulence, 6 V and 9\, would equal 45°. For 
fully anisotropic turbulence, 9 V and 9b would equal 90°, imply- 
ing that their spectra would be normal to Bo e ; . We find that, 
after t = 100, 9b is always bigger than 87°, while the velocity 
field fluctuations are found to be slightly more anisotropic than 
those in the magnetic field, 88° < 9 V < 89°, indicating that the 
turbulence in the present system is highly anisotropic, i.e., the 
turbulence is strongly confined to planes orthogonal to the guide 
magnetic field. Moreover we have verified that the dynamics is 
almost incompressible. 

The important new result of this letter refers to the thermody- 
namical behavior of the system resulting from the self-consistent 
heating due to the magnetic field-line tangling induced by pho- 
tospheric motions, once the cooling effects of conduction and 
radiation are taken into account. The fact that the total inter- 
nal energy, as mentioned above, fluctuates around a steady value 
means that the self-consistent heating mechanism due to photo- 
spheric motions is efficient enough to sustain a hot corona once 
the system has reacted to photospheric motion. Such a reaction 
leads to the formation of small scales where energy can be effi- 
ciently dissipated. The small scales, corresponding to the pres- 
ence of current sheets, are not uniformly distributed in the sys- 
tem, and as a result the temperature is spatially intermittent, as 
shown in Figure [T] where we show the midplane temperature 
contours at t — 300. 

Figure[2]shows the time evolution of T max and |j| mflA , the max- 
imum temperature and the maximum current density present in 
the loop at each time respectively. It can be seen that the T max 
temperature peaks oscillate in time and, after an initial transient, 
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Fig. 3. Temperature is organized in quasi-2D pancake-like struc- 
tures aligned with the magnetic field. The two isosurfaces for 
T = 10 6 K and T = 1.2 x 10 6 K at time t = 300 are shown re- 
spectively in yellow and grey (corresponding to T = 100, 120 in 
nondimensional units). Sample magnetic field-lines are drawn in 
red. The box has been rescaled to improve visualization. 

the T max values are well above T — 10 6 K, higher than the initial 
peak temperature of T = 8 x 10 5 K degrees. An important fact to 
emphasize is that the peaks in T max roughly coincide with bursts 
in the maximum electric current. Since currents are organized in 
sheets elongated along the mean magnetic field and the heating 
is due to the dissipation of currents we expect the heating not to 
be homogeneously distributed throughout the loop. 

This expectation is confirmed by Figure [3] where a 3-D 
snapshot of the isosurfaces of two temperatures (T = 10 6 and 
T = 1 .2 x 10 6 K) at t — 300 is presented. The spatial temperature 
morphology shows that the high temperatures are reached in se- 
lected quasi-2D regions elongated along the mean magnetic field 
where the system dynamics places at each time the current sheets 
where energy is dissipated. The higher temperatures are found 
at the center of these structures (up to 1.35 X 10 6 K at t = 300), 
while farther out the temperature falls off to a background value 
that in the central z-region averages around 6 x 10 5 K. Snapshots 
taken at different times show similar behavior with different high 
temperatures locations. 

4. Conclusions and discussion 

In this paper we have presented for the first time the thermo- 
dynamic behavior of a system where the conductive and radia- 
tive losses are balanced by the heating process due solely to 
the photospheric motions and the subsequent tangling of mag- 
netic field lines, resolving the sc ales below convection. While 
previous compressible simulation (iGudiksen & Nordlundl f2002: 
B ins ert & Peterl 1201 ll) have considered entire active regions re- 
solving at most the convective scale, it is pivotal to resolve 
smaller scales in order to understand how the plasma is heated 
and its observational consequences. 

A lot of work has been done in the past to show the reaction 
of a plasma subject to mass flow at its boundaries. It has been 



shown that the system reacts developing a weak turbulent state, 
where energy is dissipated in the current sheets resulting from 
the dynamics induced by the boundary velocity forcing. The re- 
sults of the present work show that this process leads to a loop 
where the temperature is spatially highly structured. 

A loop is therefore a multi-thermal system, where the tem- 
perature peaks around current sheets and exhibits a distribution 
of temperatures. The characteristics of this distribution (peak, 
width, background) will vary depending on the parameters of 
the loop (length, boundary velocity, mean magnetic field). 

Temperature and heating are highly nonhomogeneously dis- 
tributed so that only a fraction of the loop volume exhibits a 
significant heating at each time. Therefore emission observed in 
a spectral band origins from the small sub-volume at the cor- 
responding temperature, whose mass is considerably lower that 
the total loop mass. 

Observations currently have a spatial resolution of ~ 700 km 
and a temporal resolution of ~ 30 s. These temperature struc- 
tures are therefore unresolved, and to interpret observations cor- 
rectly the spatial and temporal intermittency must be taken into 
account. We intend to present in a following paper how our loop 
appears once the emission is averaged over the coarser observa- 
tional resolutions. At those scales the loop might result multi- 
thermal or approximately isothermal depending on the loop pa- 
rameters. 

It is very likely that the local thermodynamical equilibrium 
is lost and that the emission is due to particles accelerated in 
the current sheets rather than direct transformation of magnetic 
energy into thermal energy. This is beyond the MHD model. 

We have shown that the heating due to photospheric motions 
is able to balance the losses and that the dynamics induced by 
the such motions is a multi-scale dynamics where the energy 
is released in a number of spatial spots which account for the 
observed emission. 

The density stratification in our simulation allows limited 
evaporation from the z-boundary edges, where density is higher 
but does not reach typical chromospheric values. In future work 
we plan to imp l ement characteristic boundary conditions (e.g., 
iRappazzo et al.L 120051) that allow higher density flows through 
the z-boundaries, that can also give rise to condensation and 
prominences, and study their impact on the thermodynamics. 

Acknowledgements. This work was carried out in part at the Jet Propulsion 
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